3-freq monomials calculation

Find all 3-freq monomials up to order N, given a tolerance tol.


In [1]:
# Utilities 

from __future__ import division
import numpy as np


class MemoizeMutable:
    """
    Implement memoization with mutable arguments
    """
    def __init__(self, fn):
        self.fn = fn
        self.memo = {}
    def __call__(self, *args, **kwds):
        import cPickle
        str = cPickle.dumps(args, 1)+cPickle.dumps(kwds, 1)
        if not self.memo.has_key(str):
            self.memo[str] = self.fn(*args, **kwds)
#         else:
#             print "returning memoized calculation"

        return self.memo[str]


# got from http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays
# (faster than other alternatives)
def cartesian(arrays, out=None):
    """
    Generate a cartesian product of input arrays.

    Parameters
    ----------
    arrays : list of array-like
        1-D arrays to form the cartesian product of.
    out : ndarray
        Array to place the cartesian product in.

    Returns
    -------
    out : ndarray
        2-D array of shape (M, len(arrays)) containing cartesian products
        formed of input arrays.

    Examples
    --------
    >>> cartesian(([1, 2, 3], [4, 5], [6, 7]))
    array([[1, 4, 6],
           [1, 4, 7],
           [1, 5, 6],
           [1, 5, 7],
           [2, 4, 6],
           [2, 4, 7],
           [2, 5, 6],
           [2, 5, 7],
           [3, 4, 6],
           [3, 4, 7],
           [3, 5, 6],
           [3, 5, 7]])

    """

    arrays = [np.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = np.prod([x.size for x in arrays])
    if out is None:
        out = np.zeros([n, len(arrays)], dtype=dtype)

    m = n / arrays[0].size
    out[:,0] = np.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian(arrays[1:], out=out[0:m,1:])
        for j in xrange(1, arrays[0].size):
            out[j*m:(j+1)*m,1:] = out[0:m,1:]
    return out

The actual monomial finding functions ...


In [2]:
from collections import defaultdict, namedtuple


@MemoizeMutable
def fareySequence(N, k=1):
    """
    Generate Farey sequence of order N, less than 1/k
    """
    # assert type(N) == int, "Order (N) must be an integer"
    a, b = 0, 1
    c, d = 1, N
    seq = [(a,b)]
    while c/d <= 1/k:
        seq.append((c,d))
        tmp = int(math.floor((N+b)/d))
        a, b, c, d = c, d, tmp*c-a, tmp*d-b
    return seq


@MemoizeMutable
def resonanceSequence(N, k):
    """
    Compute resonance sequence

    Arguments:
        - N (int): Order
        - k (int): denominator of the farey frequency resonances are attached to
    """
    a, b = 0, 1
    c, d = k, N-k
    seq = [(a,b)]
    while d >= 0:
        seq.append((c,d))
        tmp = int(math.floor((N+b+a)/(d+c)))
        a, b, c, d = c, d, tmp*c-a, tmp*d-b
    return seq



def rationalApproximation(points, N, tol=1e-3, lowest_order_only=True):
    """
    Arguments:
        points: 2D (L x 2) points to approximate
        N: max order
    """
    L,_ = points.shape

    # since this solutions assumes a>0, a 'quick' hack to also obtain solutions
    # with a < 0 is to flip the dimensions of the points and explore those
    # solutions as well
    points = np.vstack((points, np.fliplr(points)))

    solutions = defaultdict(set)

    sequences = {1: set(fareySequence(1))}
    for n in range(2, N+1):
        sequences[n] = set(fareySequence(n)) - sequences[n-1]

    for h,k in fareySequence(N,1):
        if 0 in (h,k):
            continue
        # print h,k
        for x,y in resonanceSequence(N, k):

            # avoid 0-solutions
            if 0 in (x,y):
                continue

            norm = np.sqrt(x**2+y**2)

            n = np.array([ y/norm, x/norm]) * np.ones_like(points)
            n[points[:,0] < h/k, 0] *= -1  # points approaching from the left

            # nomenclature inspired in http://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation
            ap = np.array([h/k, 0]) - points
            apn = np.zeros((1,L))
            d = np.zeros_like(points)

            apn = np.sum(n*ap, 1, keepdims=True)
            d = ap - apn*n

            ## DON'T RETURN IMMEDIATELY; THERE MIGHT BE OTHER SOLUTIONS OF THE SAME ORDER OR HIGHER
            indices, = np.nonzero(np.sqrt(np.sum(d*d,1)) <= tol)
            for i in indices:
                # print "h/k:", h , "/", k
                # print "point:", points[i,:]
                if points[i,0] >= h/k:
                    if i<L:
                        # print "non-flipped >= h/k"
                        solutions[i].add((x,-y, h*x/k))
                        # print i, (x,-y, h*x/k)
                    elif x*(-y)<0:  # only consider solutions where (a,b) have different sign for the "flipped" points (the other solutions should have already been found for the non-flipped points)
                        # print "flipped >= h/k"
                        solutions[i-L].add((-y, x, h*x/k))
                        # print i-L, (-y, x, h*x/k)
                else:
                    if i<L:
                        # print "non-flipped < h/k"
                        solutions[i].add((x, y, h*x/k))
                        # print i, (x, y, h*x/k)
                    elif x*y>0:  # only consider solutions where (a,b) have different sign for the "flipped" points (the other solutions should have already been found for the non-flipped points)
                        # print "flipped < h/k"
                        solutions[i-L].add((y, x, h*x/k))
                        # print i-L, (y, x, h*x/k)

    if lowest_order_only:
        # removed = 0
        for k in solutions:
            # keep lowest order solutions only
            lowest_order = 2*N
            s = set([])
            for sol in solutions[k]:
                K = abs(sol[0])+abs(sol[1])+abs(sol[2])
                if K == lowest_order:
                    s.add(sol)
                elif K < lowest_order:
                    lowest_order = K
                    # if len(s) > 0:
                    #     print("point: ({},{}) -> removing {} for {}".format(points[k,0], points[k,1], s, sol))
                    #     removed += len(s)
                    s = set([sol])
            solutions[k] = s
        # print("Removed {} solutions".format(removed))

    return solutions


@MemoizeMutable
def threeFreqMonomials(fj, fi, allow_self_connect=True, N=5, tol=1e-10, lowest_order_only=False):
    """
    Arguments:
        fj (np.array_like): frequency vector of the source (j in the paper)
        fi (np.array_like): frequency vector of the target (i in the paper)
        N: max order
        tol: tolerance in the point to line distance calculation
    """
    from time import time
    st = time()

    # use 32bits (64 can take too much memory)
    fj = np.array([f for f in fj], dtype=np.float32)
    fi = np.array([f for f in fi], dtype=np.float32)
    Fj, Fi = len(fj), len(fi)

    cart_idx = cartesian((np.arange(Fj),
                          np.arange(Fj),
                          np.arange(Fi)))

    # we care only when y2 > y1
    cart_idx = cart_idx[cart_idx[:,1]>cart_idx[:,0]]

    if not allow_self_connect:
        cart_idx = cart_idx[(cart_idx[:,0] != cart_idx[:,2]) & (cart_idx[:,1] != cart_idx[:,2])]

    # actual frequency triplets
    cart = np.vstack((fj[cart_idx[:,0]], fj[cart_idx[:,1]], fi[cart_idx[:,2]])).T
    nr, _ = cart_idx.shape

    # sort in order to get a*x+b*y=c with 0<x,y<1
    sorted_idx = np.argsort(cart, axis=1)
    cart.sort()
    print("a) Elapsed: {} secs".format(time() - st))
    all_points = np.zeros((nr, 2), dtype=np.float32)
    all_points[:,0] = cart[:,0] / cart[:,2]
    all_points[:,1] = cart[:,1] / cart[:,2]
    # del cart
    print("b) Elapsed: {} secs".format(time() - st))

    redundancy_map = defaultdict(list)
    for i in xrange(all_points.shape[0]):
        redundancy_map[(all_points[i,0],all_points[i,1])].append(i)
    del all_points
    print("c) Elapsed: {} secs".format(time() - st))

    points = np.array([[a,b] for a,b in redundancy_map])
    print("d) Elapsed: {} secs".format(time() - st))

    exponents = rationalApproximation(points, N, tol=tol, lowest_order_only=lowest_order_only)
    print("e) Elapsed: {} secs".format(time() - st))


    monomials = [defaultdict(list) for x in fi]
    M = namedtuple('Monomials', ['indices', 'exponents'])
    # FIXME: is there a way to reduce the number of nested loops?
    for k in exponents:
        x, y = points[k,0], points[k,1]
        all_points_idx = redundancy_map[(x,y)]
        sols = exponents[k]
        for a, b, c in sols:
            for idx in all_points_idx:
                j1, j2, i = (cart_idx[idx, 0], cart_idx[idx, 1], cart_idx[idx, 2])
                reordered = (sorted_idx[idx,0], sorted_idx[idx,1], sorted_idx[idx,2])
                if reordered == (0,1,2):
                    n1, n2, d = a, b, c
                elif reordered == (0,2,1):
                    # n1, d, n2 = -a, b, c
                    n1, n2, d = -a, c, b
                elif reordered == (2,0,1):
                    # d, n1, n2 = a, -b, c
                    n1, n2, d = -b, c, a
                else:
                    raise Exception("Unimplemented order!")
                if d < 0:
                    n1, n2, d = -n1, -n2, -d
                monomials[i]['j1'] += [j1]
                monomials[i]['j2'] += [j2+len(fj)]  # add offset for fast look-up at run time (will use a flattened array)
                monomials[i]['i']  += [i+len(fj)+len(fi)]   # add offset for fast look-up at run time (will use a flattened array)
                monomials[i]['n1'] += [n1]
                monomials[i]['n2'] += [n2]
                # monomials[i]['d']  += [d]
                monomials[i]['d']  += [d-1]  # small optimization to avoid extra subtraction in every step when actually running the network

                # print i, (a,b,c), (n1, n2, d)

    print("f) Elapsed: {} secs".format(time() - st))

    # make them 2d arrays for speed up
    for i, m in enumerate(monomials):
        I = np.array([m['j1'], m['j2'], m['i']], dtype=np.int).T
        E = np.array([m['n1'], m['n2'], m['d']], dtype=np.int).T
        monomials[i] = M(indices=I, exponents=E)
    print("g) Elapsed: {} secs".format(time() - st))

    return monomials

Tests ...


In [5]:
from time import time

fc = 2.0
# num_oscs = 321
num_oscs = 11
f1 = fc * np.logspace(-2.5, 2.5, num_oscs, base=2)
f2 = fc * np.logspace(-2.5, 2.5, num_oscs, base=2)

start = time()
monomials2 = threeFreqMonomials(f1, f2, N=3, tol=1e-10)

print "Run in {:.2f} seconds".format(time()-start)

print monomials
print monomials2


a) Elapsed: 0.00114393234253 secs
b) Elapsed: 0.00126600265503 secs
c) Elapsed: 0.00227499008179 secs
d) Elapsed: 0.00242304801941 secs
e) Elapsed: 0.00374698638916 secs
f) Elapsed: 0.00499510765076 secs
g) Elapsed: 0.00541806221008 secs
Run in 0.01 seconds
[Monomials(indices=array([[ 2, 15, 22],
       [ 0, 13, 22]]), exponents=array([[-1,  1,  1],
       [-1,  1,  0]])), Monomials(indices=array([[ 3, 16, 23],
       [ 1, 14, 23]]), exponents=array([[-1,  1,  1],
       [-1,  1,  0]])), Monomials(indices=array([[ 0, 13, 24],
       [ 0, 15, 24],
       [ 4, 17, 24],
       [ 2, 15, 24]]), exponents=array([[ 2,  1,  1],
       [-2,  1,  0],
       [-1,  1,  1],
       [-1,  1,  0]])), Monomials(indices=array([[ 1, 14, 25],
       [ 1, 16, 25],
       [ 5, 18, 25],
       [ 3, 16, 25]]), exponents=array([[ 2,  1,  1],
       [-2,  1,  0],
       [-1,  1,  1],
       [-1,  1,  0]])), Monomials(indices=array([[ 2, 15, 26],
       [ 0, 13, 26],
       [ 2, 17, 26],
       [ 6, 19, 26],
       [ 4, 17, 26]]), exponents=array([[ 2,  1,  1],
       [ 2,  1,  0],
       [-2,  1,  0],
       [-1,  1,  1],
       [-1,  1,  0]])), Monomials(indices=array([[ 3, 16, 27],
       [ 1, 14, 27],
       [ 3, 18, 27],
       [ 7, 20, 27],
       [ 5, 18, 27]]), exponents=array([[ 2,  1,  1],
       [ 2,  1,  0],
       [-2,  1,  0],
       [-1,  1,  1],
       [-1,  1,  0]])), Monomials(indices=array([[ 4, 17, 28],
       [ 2, 15, 28],
       [ 4, 19, 28],
       [ 8, 21, 28],
       [ 6, 19, 28]]), exponents=array([[ 2,  1,  1],
       [ 2,  1,  0],
       [-2,  1,  0],
       [-1,  1,  1],
       [-1,  1,  0]])), Monomials(indices=array([[ 5, 18, 29],
       [ 3, 16, 29],
       [ 5, 20, 29],
       [ 7, 20, 29]]), exponents=array([[ 2,  1,  1],
       [ 2,  1,  0],
       [-2,  1,  0],
       [-1,  1,  0]])), Monomials(indices=array([[ 6, 19, 30],
       [ 4, 17, 30],
       [ 6, 21, 30],
       [ 8, 21, 30]]), exponents=array([[ 2,  1,  1],
       [ 2,  1,  0],
       [-2,  1,  0],
       [-1,  1,  0]])), Monomials(indices=array([[ 7, 20, 31],
       [ 5, 18, 31]]), exponents=array([[2, 1, 1],
       [2, 1, 0]])), Monomials(indices=array([[ 8, 21, 32],
       [ 6, 19, 32]]), exponents=array([[2, 1, 1],
       [2, 1, 0]]))]

In [ ]: